Binary black hole mergers in magnetized disks: simulations in full general relativity 
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We present results from the first fully general relativistic, magnetohydrodynamic (GRMHD) sim- 
ulations of an equal-mass black hole binary (BHBH) in a magnetized, circumbinary accretion disk. 
We simulate both the pre and post-decoupling phases of a BHBH-disk system and both "cooling" 
and "no-cooling" gas flows. Prior to decoupling, the competition between the binary tidal torques 
and the effective viscous torques due to MHD turbulence depletes the disk interior to the binary 
orbit. However, it also induces a two- stream accretion flow and mildly relativistic polar outflows 
from the BHs. Following decoupling, the accretion rate is reduced, while the EM luminosities peak 
near merger due to shock heating. This investigation, though preliminary, previews more detailed 
GRMHD simulations we plan to perform in anticipation of future, simultaneous detections of grav- 
itational and electromagnetic radiation from a merging BHBH-disk system. 
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When galaxies merge, massive black hole binaries 
(BHBHs) are likely to form in gas rich environments pQ. 
These systems may provide a unique opportunity for a si- 
multaneous detection of electromagnetic (EM) and grav- 
itational waves (GWs) from the same cosmic event. 

If BHBHs are embedded in a gas with negligible an- 
gular momentum, the accretion flow will resemble the 
Bondi or Bondi-Hoyle-Lyttleton accretion solutions [2]- 
H]. When the gas has intrinsic angular momentum, it 
will form a circumbinary disk that accretes via angular 
momentum transport induced by an effective viscosity. 
The BHBH-disk problem has been studied extensively 
in Newtonian gravitation, both analytically [SH9], where 
the emphasis has been on low mass-ratio binaries, and 
numerically [5] [TQl [TT] where equal-mass systems have 
been the focus. The goal is to compute observable EM 
"precursor" and "aftermath" radiation that will accom- 
pany the GW signal. While the vacuum BHBH [12] and 
accretion onto a single BH problems in GR are now well 
developed, simulations of disk accretion onto BHBHs are 
still in their infancy [2H4} (13] [T3], Vacuum calculations 
offer an accurate description of the spacetime and the 
GWs emitted close to merger in typical cases. Determin- 
ing the EM signatures requires a GRMHD computation 
in this dynamical BHBH spacetime. 

Here we report the first fully GRMHD simulation of 
a magnetized, circumbinary BHBH accretion disk. The 
effective viscosity driving accretion arises from MHD tur- 
bulence triggered by the magnetorotational instability 
(MRI) [15]. This effective viscosity competes with the 
tidal torques exerted by the binary, so that a quasi- 
stationary state is reached prior to binary-disk decou- 
pling [6j [7] . This state has been simulated both in New- 
tonian [6 and in Post-Newtonian [16] gravitation. Typi- 



* bfarris2@illinois.edu 

t Also at Department of Astronomy & NCSA, University of Illinois 
at Urbana-Champaign, Urbana, IL 61801 



cally, the computational domain excludes the region near 
the BHs and artificial inner boundary conditions are im- 
posed. Recent Newtonian studies [17] make clear the 
importance of imposing the correct boundary conditions 
on the flow inside the central disk hollow and near the 
BHs, further motivating a treatment in full, dynamical 
GR whereby BH horizons can be modeled reliably. 

The basic evolution of the system is as follows: For 
large binary separations a, the inspiral time due to GW 
emission is much longer than the viscous time (£gw ^ 
t v i S ), so that the disk settles into a quasi-stationary state. 
For equal-mass BHs, the binary tidal torques carve out 
a partial hollow in the disk [5j [6] [8] [10] of radius ~ 2a 
and excite spiral density waves throughout the disk, that 
dissipate and heat the gas. However, gas can pene- 
trate the hollow in response to the time-varying tidal 
torque [8j [lOj [16] [18] . At sufficiently small separations 
tew ^ ^vis? and the BHBH decouples from the disk. 
The disk structure at decoupling crucially determines its 
subsequent evolution and the EM emission. GW emis- 
sion close to merger leads to mass loss [19] [20] and may 
induce remnant BH recoil [21], which give rise to fur- 
ther characteristic EM signatures. Here we simulate the 
system in two different epochs: (I) The pre-decoupling 
phase (tow > Uis) and (II) the post-decoupling phase 
(tew < tvis), including the inspiral and merger of the 
BHBH. We consider equal-mass, nonspinning binaries. 
While the BH mass scales out, we are primarily inter- 
ested in total (ADM) masses M > 1O 6 M and low den- 
sity disks for which the tidally-induced binary inspiral 
and the disk self- gravity are negligible. 

We use the Illinois numerical relativity code to carry 
out our simulations. The code has been extensively 
tested [22] [23] and used in our earlier BHBH simula- 
tions in gaseous media [2] [18]. For details and equa- 
tions see [22H24] . The main new feature concerns our 
vector potential (A^ = Qn^ + A^) formulation for the 
magnetic induction equation, where n M is the future- 
directed timelike unit vector normal to a t =const. slice 
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and n^A^ = 0. We introduce a new generalized Lorenz 
gauge condition V M ^4 M = (.n^A^, where £ is a parameter 
(typically £ = 4/M) and M is the total BHBH (ADM) 
mass. This condition results in damped traveling EM 
gauge modes, preventing the spurious emergence of IB- 
fields associated with refinement boundaries [25] . 

The disk initial data represent an equilibrium disk or- 
biting a single Schwarzschild BH [18] [26] with an inner 
disk edge at R[ n = 18M = 1.8a, where the specific an- 
gular momentum i- 1Ti = 5.15M at R ln , and a nearly Ke- 
plerian rotation profile parameter q = 1.7. We adopt a 
T-law equation of state with T = 5/3, appropriate for 
a disk composed of an ideal, nonrelativistic gas. We 
seed the disk with a weak poloidal B-field as described 
in [24]. The maximum relative strength of the initial IB- 
field in the equatorial plane is (Pm / P)max — 0.025. Here 
Pm = B 2 /8tt is magnetic pressure, P is gas pressure, and 
B^ is the magnetic field measured in the comoving frame 
of the fluid. 

Prior to decoupling we can neglect the slow BHBH 
inspiral. We model the spacetime during this epoch by 
adopting the BHBH metric derived in the conformal thin- 
sandwich (CTS) formalism [27], whereby the spacetime 
is stationary in the corotating frame (see [18] for details). 
The inner part of the disk settles into a quasiequilibrium 
state on a "viscous" time scale 

where v is the effective viscosity induced by MHD tur- 
bulence [15]. This viscosity can be fit (approximately) 
to an 'a-disk' law for purposes of analytic estimates. 
Here R is the disk radius, v{R) = (2/3)a(P/ Po)^ 1 « 
(2/3)a(R/M)^ 2 (H/R) 2 M, H is the disk scale height, 
and we have assumed vertical hydrostatic equilibrium 
to derive an approximate relationship between Pj po and 
H/R (see [28]). Equating the viscous time scale and the 
GW inspiral time scale yields the decoupling separation 

where the normalizations give the typical parameters our 
simulations obtain for the relaxed state. Note, that for 
the geometrically thick, magnetic disks we treat, the ex- 
pected decoupling radius is an order of magnitude smaller 
than typical thin-disk cases [6] [29]. We thus set our 
initial binary separation at a/M = 10 (orbital period 
2n/n = 225M). 

We evolve the system using the CTS spacetime for 
~ 45 binary orbits (10,000 M) to allow the inner parts 
of the disk to settle into a quasistationary state. This 
epoch (1) models the pre-decoupling phase and (2) pro- 
vides realistic, relaxed disk initial data for the post- 
decoupling inspiral phase. We model the post-decoupling 
phase by continuing the GRMHD evolution in the dy- 
namical spacetime of the inspiraling and merging BHBH 



12 3 4 




R / a 

FIG. 1. The disk a parameter (upper panel) and surface- 
density £ (lower panel) profiles. So is the maximum surface 
density at t = 0. Dotted magenta line is the initial data, solid 
lines are at decoupling, and dashed lines are at merger. Black 
lines are from the "no-cooling" case, and blue lines are from 
the "cooling" case. 



binary. We treat two extreme opposite limiting cases: 
"no-cooling" , which allows for gas heating via shocks in- 
duced by tidal torques and MHD turbulence, and "cool- 
ing" , which removes all the heat generated via an effective 
local emissivity A of the form T Mzy ;zy = —h.u^ as in [30] . 
Our "cooling" case, though artificial, provides a repre- 
sentative example of the effects of cooling and has been 
adopted in previous work (e.g. [I6j[3l]). The particular 
"cooling" prescription we use drives the gas to isentropic 
behavior, i.e. P/Pq = const. The cooling timescale is 
set to the local, Keplerian orbital period. Our simula- 
tions resolve the BH horizons and we impose no inner 
boundary conditions. In the pre-decoupling phase our 
grid consists of a hierarchy of 6 refinement levels with 
(coarsest, finest) resolution of (5.33M, 0.16M) and outer 
boundary at 250M. We resolve the fastest-growing MRI 
mode by at least 10 grid points in the bulk of the inner 
disk. We add two extra levels centered on each BH in the 
post-decoupling phase, increasing the (coarsest, finest) 
resolution to (4M, M/32). After merger and ringdown 
we freeze the spacetime evolution, but continue to evolve 
the plasma. Equatorial symmetry is imposed through- 
out. We normalize results to those for a single BH that 
we evolved with the same initial magnetized disk and BH 
mass equal to M. 

The initial disk (see Fig. [I]) is not in equilibrium around 
the binary as it is perturbed by the binary torques. This 
leads to spiral density waves in the disk which dissi- 
pate and heat the gas. The gas gains angular momen- 
tum and the surface density profile moves slightly out- 
ward. Magnetic winding converts the poloidal field into 
one with a large toroidal component. MRI is induced, 
resulting in turbulent flow. After about 20 binary or- 
bits (~ 4 — 5 disk orbits at the pressure maximum) the 
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FIG. 2. Orbital plane snapshot of rest-mass density 
log (po/po, max) from the "no-cooling" simulation at t ~ 
10000M in the relaxed disk, prior to decoupling. The inset 
zooms in on the region close to the BHs. 
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FIG. 4. Meridional snapshot of magnetic pressure 

log(PM / po,max) and fluid velocity vectors at t ~ 10000M for 
the "no-cooling" case. 
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FIG. 3. Time-averaged binary accretion rate Mbhbh, nor- 
malized to the average value for a single BH Mbh, versus 
time. Colors have the same meaning as in Fig. [I] Left 
panel: pre-decoupling (a =const) phase. Right panel: post- 
decoupling (inspiral) phase. Mbh = O.45ni2M|M0yr _1 
where M 8 M/10 8 M o and n 12 n/10 12 cm~ 3 is the ini- 
tial maximum gas particle number density. 



MRI saturates, driving disk accretion onto the BHBH. In 
the relaxed disk prior to decoupling we measure a time- 
averaged Maxwell-stress as in [32 at 20M < R < 30M, 
and find a = 0.13 for the "no-cooling" (see Fig. [T]) and 
a = 0.2 for the "cooling" case. The magnetic-to-gas- 
pressure ratio ranges from 0.1 to 5 in the bulk of the 
disk. 

Cooling influences the global disk structure. In partic- 
ular, we observe matter pile- up near the inner disk edge 
only with cooling (see Fig. [I]), as has previously been 
found in [7J [TO] [TT] [16] . The binary maintains a partial 
hollow in the disk (see Fig. [T]) by exerting torques on the 
plasma, while the MRI-induced effective viscosity drives 
matter inward. Cooling leads to smaller scale height and 
lower which explains the enhanced pile- up at Ri n . We 
confirm the result in [TO] [TT] [TO] [18] that non-negligible 
amounts of gas are present inside the cavity. 

Accretion occurs predominantly via two spiral density 



streams inside the cavity (see Fig.j2|. We find that accre- 
tion exhibits an alternating pattern by accreting primar- 
ily on one of the BHs for about half a binary orbit. This is 
similar to flow features observed in [11] . The behavior has 
been attributed to a gradual increase of disk-eccentricity 
[TO] , which weakens one of the two streams when the BHs 
pass near the disk-apocenter and strengthens the other 
stream at pericenter. 

Prior to decoupling, but after an initial transient phase 
(5000M <t< tvis(^in))? the accretion rate Mbhbh set- 
tles to values comparable to those onto a single BH of 
mass M (see Fig. [3| . We perform a Fourier analysis of 
Mbhbh and find that the strongest contributions arise 
at / ~ 2/3Q in both of our cases. This is likely asso- 
ciated with the dominant (2,3) Lindblad resonance [10 . 
During the early inspiral the inward drift of the disk edge 
lags behind the binary orbital decay, decreasing Mbhbh- 
In contrast to the magnetic-free case [18], the accretion 
streams remain present until merger at t m for the "no- 
cooling" case and up until t — t m > — 400M for the 
"cooling" case. At merger Mbhbh decreases gradually 
to about 30% of the single, quasi-stationary BH accre- 
tion rate, Mbh (see Fig. 31). 

The remnant BH settles down via quasi-normal mode 
ringing to a Kerr-like BH with mass Mf ~ 0.95M and 
dimensionless spin s = Jf/Mj = 0.68. In the "cool- 
ing" case the inner disk edge reaches Ri n (tm) °° 10M 
at merger, while the edge is more dispersed in the "no- 
cooling" case (see Fig. [I]). 

Prior to decoupling we detect persistent, magnetized, 
mildly relativistic (v > 0.5c) collimated outflows in the 
polar regions (see Fig. |4| in both cases. After merger 
there is an increase of the velocities in the outflow reach- 
ing Lorentz factors W < 4 for the "no-cooling" case, 
while the outflow in the "cooling" case remains rela- 
tively unchanged. These properties persist throughout 
the postmerger evolution. After merger the effective, tur- 
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FIG. 5. The Poynting luminosity Lem measured at R = 
100M as a function of time (blue line for "cooling", black 
line for "no-cooling"). The cooling luminosity L coo i from the 



"cooling" case (red line). Mbhc 2 = 2.58 • 10 46 ergs" 
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bulent viscous torque will cause the gas to refill the cavity 
and accrete on the merger remnant. Thus, brightening 
crucially depends on the surface density at decoupling. 
If there is only a small pile-up and the majority of gas 
lies at radii R > 40M, then a significant brightening 
will take t vi8 (R > 40M) > (9(10 4 M) following merger. 
Using Eq. Q we estimate t vis (R = 40M) > 22000M 
(a ~ 0.13, H/R ~ 0.3) for the "no-cooling" case and 
t vis (R = 14M) > 6800M (a ~ 0.27, H/R = 0.17) for 
the "cooling" case, respectively. We observe the sur- 
face density profile diffusing inward, but we have not 
followed the evolution for this long. The Poynting lu- 
minosity (Lem = j —T r [ EM ^ y/^gdS) across a sphere at 
R = 100M and the (approximate) luminosity from disk 
cooling (L coo i = j kutyf^dPx) are shown as functions 
of time in Fig. [5j We find a sudden enhancement in the 
Poynting luminosity for the "no-cooling" case beginning 
at merger and growing for < 200M. This enhancement 
originates from shocked gas in the immediate vicinity of 
the binary and escapes mainly through the polar regions. 
It is absent for the "cooling" case, for which there is less 
gas near the binary prior to merger. Heat generated by 
shocks is also removed in the cooling case, leaving less 
energy available for conversion to a Poynting flux. In 
Fig. ([5| we also plot the cooling luminosity for the "cool- 
ing" case, for which we also observe a sudden jump at 
merger. We find that the outward flux of kinetic en- 
ergy is much smaller in both cases. Total efficiencies 
e = L/Mbhc 2 increase to e = 0.25 at merger. 

The Poynting luminosity is presumably reprocessed at 
larger distance from the remnant. The characteristic fre- 
quencies of the total emitted EM radiation will depend 
on the BH masses, disk densities and dominant cooling 
mechanisms. We plan to investigate these issues, along 
with other precursor (e.g. twin jets) and afterglow ef- 
fects, different mass ratios, and different BH spins more 
thoroughly in future work. 
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